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Kent,  Washington 


A  stress  analysis  method  is  presented  which  includes  the  capa¬ 
bility  for  transient,  non  homogeneous  temperature  distribution.  The 
method  is  based  on  the  following  assumptions:  linear  viscoelasticity 
with  hereditary  integral  form  of  stress- strain  relation;  validity  of  the 
reduced  time  hypothesis:  bulk  modulus  constant  in  time;  homogeneous, 
isotropic  material.  The  associated,  finite  element  computer  program, 
called  TAP  (Thermoviscoelastic  Analysis  Program),  is  aplane  strain 
formulation  with  uniform  strain  and  linear  temperature  distribution 
assumed  in  each  element.  The  element  matrices  involve  superposition 
integrals  which  are  approximated  numerically  for  stepping  out  the 
time  varying  solution.  TAP  solutions  for  simple  problems  are  com¬ 
pared  with  exact  and  other  approximate  solutions  and  with  other 
approximate  methods  to  more  complex  problems.  Different  numerical 
procedures  are  considered,  in  particular  one  that  avoids  integration 
back  to  the  time  origin  thus  minimi  zing  computer  storage  and  solution 
times. 


*  Research  Engineer,  Missile  and  Information  Systems  Division,  Structural  Analysis  Group. 
The  assistance  of  Messrs.  J.  Murray,  J.  Loski,  and  M.  Wolff  of  Structures  R&D  Computing 
is  acknowledged  in  developing  computer  codes. 


489 


Report  Documentation  Page 

Form  Approved 

OMB  No.  0704-0188 

Public  reporting  burden  for  the  collection  of  information  is  estimated  to  average  1  hour  per  response,  including  the  time  for  reviewing  instructions,  searching  existing  data  sources,  gathering  and 
maintaining  the  data  needed,  and  completing  and  reviewing  the  collection  of  information.  Send  comments  regarding  this  burden  estimate  or  any  other  aspect  of  this  collection  of  information, 
including  suggestions  for  reducing  this  burden,  to  Washington  Headquarters  Services,  Directorate  for  Information  Operations  and  Reports,  1215  Jefferson  Davis  Highway,  Suite  1204,  Arlington 

VA  22202-4302.  Respondents  should  be  aware  that  notwithstanding  any  other  provision  of  law,  no  person  shall  be  subject  to  a  penalty  for  failing  to  comply  with  a  collection  of  information  if  it 
does  not  display  a  currently  valid  OMB  control  number. 

1.  REPORT  DATE 

OCT  1968 

2.  REPORT  TYPE 

3.  DATES  COVERED 

00-00-1968  to  00-00-1968 

4.  TITLE  AND  SUBTITLE 

Finite  Elements  in  Linear  Viscoelasticity 

5a.  CONTRACT  NUMBER 

5b.  GRANT  NUMBER 

5c.  PROGRAM  ELEMENT  NUMBER 

6.  AUTHOR(S) 

5d.  PROIECT  NUMBER 

5e.  TASK  NUMBER 

5f.  WORK  UNIT  NUMBER 

7.  PERFORMING  ORGANIZATION  NAME(S)  AND  ADDRESS(ES) 

Air  Force  Flight  Dynamics  Laboratory, Wright  Patterson  AFB, OH, 45433 

8.  PERFORMING  ORGANIZATION 

REPORT  NUMBER 

9.  SPONSORING/MONITORING  AGENCY  NAME(S)  AND  ADDRESS(ES) 

10.  SPONSOR/MONITOR'S  ACRONYM(S) 

11.  SPONSOR/MONITOR'S  REPORT 
NUMBER(S) 

12.  DISTRIBUTION/AVAILABILITY  STATEMENT 

Approved  for  public  release;  distribution  unlimited 

13.  SUPPLEMENTARY  NOTES 

See  also  AD0703685,  Proceedings  of  the  Conference  on  Matrix  Methods  in  Structural  Mechanics  (2nd) 

Held  at  Wright-Patterson  Air  Force  Base,  Ohio,  on  15-17  October  1968. 

14.  ABSTRACT 

15.  SUBIECT  TERMS 

16.  SECURITY  CLASSIFICATION  OF: 

17.  LIMITATION  OF 
ABSTRACT 

18.  NUMBER 
OF  PAGES 

28 

19a.  NAME  OF 
RESPONSIBLE  PERSON 

a.  REPORT 

unclassified 

b.  ABSTRACT 

unclassified 

c.  THIS  PAGE 

unclassified 

Standard  Form  298  (Rev.  8-98) 

Prescribed  by  ANSI  Std  Z39-18 


AFFDL-TR-6  8-150 


SECTION  I 
INTRODUCTION 

This  paper  was  motivated  by  problems  of  solid  propellant  grain  stress  analysis  that  arise 
under  transient  thermal  environment.  At  present,  most  stress  analyses  of  rocket  grains  use  an 
elastic  analysis  which  by  various  techniques  is  made  to  approximate  the  viscoelastic  character 
of  the  solid  propellant  response. 

One  such  technique  uses  elastic  moduli  equal  to  the  corresponding  viscoelastic  properties 
at  the  desired  time  and  temperature.  This  method,  commonly  referred  to  as  quasi-elastic, 
essentially  ignores  the  entire  past  history  of  loading  and  environment  and  may  be  a  gross 
approximation  to  the  true  response. 

A  generally  more  preferable  technique  is  based  on  the  correspondence  principle,  or 
elastic- viscoelastic  analogy,  first  elucidated  by  Alfrey  (Reference  1)  which  permits  one,  using 
an  elastic  solution,  to  obtain  the  corresponding  viscoelastic  solution  to  the  same  geometric 
problem  with  equivalent  time-step  boundary  conditions.  The  procedure  is  exact  if  a  closed 
form  elastic  solution  exists  and  if  the  inverse  of  the  Laplace  transformed  viscoelastic  solution 
can  be  obtained  exactly.  For  use  with  numerical  solutions,  Schapery,  (Reference  2)  has  shown 
how  the  method  is  applied  as  an  approximate  inverse  Laplace  transform  technique.  The 
technique  is  appropriate  where  the  required  correspondence  exists  between  the  elastic  and 
transformed  viscoelastic  problem. 

As  pointed  out  by  others,  Morland  andLee  (Reference  3),  the  correspondence  principle  fails 
to  hold  for  nonhomogeneous,  transient  temperature  distribution,  although  still  valid  for  either 
homogeneous,  transient  temperature  or  nonhomogeneous,  steady  state  temperatures.  The 
difficulty  involves  one  of  two  things,  depending  on  whether  the  problem  is  formulated  in  real 
time  or  shifted  time.  Interms  of  real  time,  the  stress-strain  law  is  not  a  convolution  integral; 
if  the  shifted  time  is  used,  the  transformed  viscoelastic  equilibrium  and  strain  displacement 
equations  are  not  equivalent  to  the  corresponding  elastic  equations. 

To  circumvent  these  problems,  Hilton  and  Russell  (Reference  4)  impose  conditions  of 
constant  temperature  over  time  increments  and  apply  the  correspondence  principle  on  an 
incremental  basis.  This  method  has  been  used  successfully  by  Valanis  and  Lianis  (Reference  5) 
although  its  application  becomes  somewhat  involved. 
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Regardless  of  what  indirect  means  are  used  to  obtain  a  solution  to  the  thermoviscoelastic 
problem,  one  is  still  faced  with  assessing  the  accuracy  of  the  solution,  given  a  particular 
viscoelastic  material,  history  of  load  and  environment,  and  solution  time  realm.  Therefore  a 
need  exists  for  a  direct  method  of  attack  on  the  problem  -  one  in  which  there  is  confidence  of 
an  accurate  solution  if  a  fine  enough  time  grid  is  used  and  general  enough  to  include  two 
dimensional,  irregular  geometries  and  arbitrary  material  properties. 

A  first  step  in  this  direction  was  taken  by  Taylor  and  Chang  (Reference  6)  who  formulated 
the  one  dimensional  plane  strain  circular  cylinder  problem  by  a  finite  element  procedure, 
then  demonstrated  its  application  to  steady  state  temperature  fields.  As  they  point  out,  com¬ 
puter  time  can  quickly  become  excessive  for  large  problems,  which  is  due  partly  by  the  need 
for  excessive  out-of-core  storage.  The  solution  procedure  involves  generation  and  solution  of 
a  new  set  of  stiffness  equations  for  each  time  step.  The  solutions  are  saved  for  all  time  steps 
and  used  to  generate  new  load  columns  which  involve  summation  back  to  the  time  origin,  thus 
requiring  a  great  amount  of  storage  for  long-time  realms  and  large  numbers  of  freedoms. 

The  present  paper  is  a  two  dimensional,  plane  strain,  finite  element  formulation  of  the 
thermoviscoelastic  problem,  in  which  steps  are  taken  to  overcome  some  of  the  aforementioned 
difficulties.  It  is  a  direct  approach,  not  relying  on  equivalent  elastic  solutions  or  transform 
methods. 

The  assumption  of  time-constant  bulk  modulus  is  made  (an  apparent  good  assumption  for 
some  solid  propellant  materials  for  example,  Reference  7)  which  leads  to  a  particularly  simple 
time-varying  form  of  the  stiffness  matrix  and  minimizes  regeneration  of  matrices  at  each 
time  step.  Also,  it  is  shown  how  the  history  or  memory  effects  can  be  formulated  so  that  the 
memory  load  (a  generalized  load  which  accounts  for  the  history  of  deformation)  may  be  ex¬ 
pressed  as  a  recurrence  relation  from  one  time  step  to  the  next  rather  than  requiring  sum¬ 
mation  to  the  time  origin. 
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SECTION  n 

ONE  DIMENSIONAL  PROTOTYPE  EQUATION 


Certain  aspects  of  the  general  equations  in  the  next  section  can  be  illustrated  in  one  dimen¬ 
sion.  We  therefore  consider  briefly  the  numerical  inversion  of  the  uniaxial,  viscoelastic 
stress-strain  relation. 


cr{  t ) 


£10 ),(,)-/  .<>')  d.' 

0* 


dt  * 


(I) 


where  cr  (t)  and  e  (t)  are  uniaxial  stress  and  strain.  E(t)  is  the  uniaxial  stress  relaxation 
modulus  defined  as  the  stress  response  to  a  unit  step  strain.  In  Equation  1  explicit  allowance 
for  an  initial  step  strain  has  been  made.  When  cr  (t)  is  prescribed  and  E(t)  is  given  as  a 
tabulated  numerical  function,  it  is  desired  to  find  a  numerical  solution,  e  (t).  When  cr(t)  is 
a  unit  step,  the  solution  e  (t)  is  the  uniaxial  creep  compliance,  D(t),  The  numerical  form  of 
Equation  1  using  a  trapezoidal  approximation  is 


i=k-'  i 

cr(.k)=E(0lc(tk)  -  X  - - - [E  >-E( Vi’] 

i  =  I 


(2) 


or  rearranging, 

rEtO)+E(t1,-lk_|) 


1  €  U,  -  |  ^  r  l 

<7- »„)+  - —  [E(0)-E(tk-tk.,)J 


+  '  X  «*U|  I  [El'k-'i  +  l  !]  =  O'  I  tk»  +V(tk) 

i  -  I 


(3) 


K  =  I,  2  ,3, 


t,  =  0 


e(t .  .  , )  +  e (t.  ) 

,  i  +  l  1 

e  (t.  )  -  - r 


from  which  is  determined  successively  €  (t^).  The  “memory”  term,  V ( t^) ,  which  involves 
the  past  solutions,  is  known.  Lee  and  Rogers  (Reference  8)  have  discussed  convergence  using 
the  trapezoidal  approximation. 


A  simpler  approximation  is  to  assume  E(t)  constant  over  each  time  step,  a  staircase 
function  falling  below  the  exact  E{t).  The  advantage  of  this  approximation  is  that  an  upper  value 
for  the  time  varying  solution  (creep  compliance)  is  obtained,  which  can  be  demonstrated 
analytically.  The  trapezoidal  approximation  on  the  other  hand  appears  to  give  a  value  below 
the  exact  solution.  This  information  could  possibly  be  used  to  bracket  the  solution  to  certain 
types  of  problems. 
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In  order  to  gain  some  insight  into  the  numerical  requirements  for  differing  viscoelastic 
materials,  two  representative  relaxation  functions  were  inverted,  one  a  very  simple,  standard 
linear  solid  representation  and  the  other  a  more  complex,  broad  band  representation.  The 
results  are  shown  in  Figures  1  and  2.  The  inversion  for  the  simpler  material  converged 
rapidly  but  the  reflected  curve  1/E  (t)  deviated  greatly  from  D(t).  The  more  complex  represen¬ 
tation,  for  which  a  closed  form  exact  creep  compliance  is  unobtainable,  showed  much  slower 
convergence  but  an  apparently  reasonable  representation  of  D(t)  by  1/E  (t).  We  might  therefore 
expect  quasi-elastic  solutions  to  be  acceptable  engineering  approximation  in  some  cases. 


Ordinarily  the  summation  in  Equation  2  must  be  carried  out  back  to  the  time  origin  re¬ 
quiring  that  every  previous  solution  be  saved.  In  the  multi-dimensional  case  considered  in  the 
next  section,  vector  solutions  must  be  saved  and  summed  at  each  time  step  so  that  storage 
requirements  and  perhaps  computer  running  time  can  become  critical  for  large  t.  There  are 
ways  to  avoid  this  problem,  none  of  which  in  general  appear  to  be  completely  satisfactory. 


One  way  is  to  truncate  the  summation  fori  <<  k  when  the  terms  j^E(tk~t.)  -  Eft^-t,  2)j 
€  (t.  i)  become  negligibly  small  since  E(t)  is  a  monotonically  decreasing  function.  In  the 
absence  of  temperature  effects  [  E(tk-t.)  “  E^k-ti-2^]  Senerally  decreases  rapidly;  how¬ 
ever,  a  past  history  of  relatively  large  deformation  could  offset  this  so  that  care  must  be 
exercised.  Further,  the  effect  of  low  temperature  is  to  retard  the  rate  of  decrease  of  E(t). 


Another  way  of  avoiding  the  storage  requirement  is  to  replace  the  summation  by  a  re¬ 
currence  relation.  Zak  (Reference  9)  has  shown  how  this  can  be  done  by  expressing  the  relaxa¬ 
tion  modulus  as  a  Prony  series  as  follows; 


E  ( 


>  *A 


The  Expression  2  then  becomes 


j  =q  — f~ 

I  Aj  e  J 


(4) 


cr  (tR  )  =  E(0)e(tk  )  - 


k-2 
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CREEP  COMPLIANCE  ,  D  (t) 
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Figure  1.  Numerical  Inversion  of  the  Relaxation  Modulus 
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Interchanging  the  summation,  the  last  term  is  written  as 

q 


k-2  _lk  JlL  q 

Z  A  I  e  Ti  (  «  r/  -e  rj  )  €*(t:J  =  I  A .  a  •  , 


J-.  J  i=. 


i=  i 


i  j  *  k 


(6) 


where 
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a  *  e  i  (e  J  -e  i  )« *"  (t  He  J 
j.k  k-2 

li+l  ii_ 

(e  rj  -e  rj  )  e*(tj  ) 


- 1 

k-3 


k-i 


1  •  ri 

i  =  i 


(7) 


or 


-ttk-tk-.  > 

r. 


a  =  e 

I,  k 


-(t.  -  t  ) 

k-i  k-2 

rj 


j ,  k  -i  ] 


(8) 


The  solution  of  Equation  1  may  now  be  written 


fttk’:  E(0)  +  E(tk  -tk.i  ) 


->  [CTl,k l+ 


El0>- 


-  ‘V, )+^  Ai  <vl 


(9) 


By  using  Zak’ s  method  we  have  exchanged  the  summation  over  the  time  realm  for  a  sum¬ 
mation  over  the  Prony  series  terms  plus  a  recurrence  relation.  Only  the  immediately  past 
two  solutions  are  now  required  to  be  saved  to  account  for  memory  effects.  The  adequacy  of 
Zak* s  method  depends  on  how  accurately  the  Prony  series  represents  the  relaxation  modulus. 
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SECTION  HI 

FINITE  ELEMENT  DERIVATIONS 

STRESS-STRAIN  RELATIONS 


The  thermoviscoelasticity  theory  used  in  the  finite  element  derivations  has  been  pre¬ 
sented  elsewhere  (Reference  10)  and  only  the  essentials  are  given  here. 


Assuming  constant  bulk  modulus,  K,  and  constant  coefficient  of  thermal  expansion,  a  ,  the 
thermoviscoelastic  constitutive  relations  may  be  written,  using  cartesian  tensor  notation,  as 


(  X,t  )  = 


2G  ( o) 

|J 


(  X  ,t  )  -  2 


l 


ao_<€-e'> 

at 


ij 


(  X  /)dt' 


+  8..  r  K  -  -|  G(o  >1  €  (X,f)+S  —  f  ]  (X  ,/)  ' 

'J  L  3  J  femm  ’  °ij  3  J  St  €mm  ’  'dt 

o+ 

-8  ..  3a  K  [t  (X,t  )  -  T  (  X  ,  o )] 


(10) 


where  allowance  for  an  initial  response  has  been  made  explicit,  G(t)  is  the  shear  relaxation 
modulus  given  for  an  arbitrary  “material  reference  temperature”.  To,  The  notation  X 
denotes  spatial  dependence.  The  shifted  time  £  is  related  to  real  time  t  through  the  relation 


t  -€(x,t)  =  f 

o 


dt' 

at  [nx.t')] 


(ID 


where  T(X  .t)  is  temperature,  while  T^(  X  ,o)  =  structural  reference  temperature,  an  initial 
steady- state  temperature  distribution  with  associated  time  constant  stress  and  displacement 
states.  The  form  to  be  used  here  for  the  shift  function  AT(T)  makes  use  of  the  so-called  WLF 
equation  (Reference  11) 


or 


log(0  AT(T) 


C,  (T  -  To  ) 
C2  -  (  T-  To) 


=  - h  (  T) 


(12) 


(13) 
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where  C  and  C  are  constants  for  a  particular  material.  Upon  the  restriction  to  plane  strain, 
X  2 

Equation  10  becomes  in  matrix  notation 


<T(X,t)=  L  «(X,t)  +C  €(X,t)  -  3Ka  I  ®(X,t) 


where 


L  «  <X,t)  J 


dG  (£-£') 

d  r 


-4  2 


*x  (X,t’) 
€y  (X,t') 


K  — G  (  o  ) 
K  +  -j  G  (o  ) 
0 


®(X,0=  [  T  t  X  ,t)  -T  (X  ,o)] 


(14) 


(15) 


One  should  be  reminded  thatthe  law  Equation  10  is  nonlinear  in  temperature,  thus  thermal 
stress  solutions  cannot  be  superimposed.  Since  temperatures  are  prescribed  this  nonlinearity 
does  not  create  any  special  problem. 


ELEMENTAL  EQUATIONS 

Derivation  of  the  elemental  equations  proceeds  from  the  principle  of  virtual  displace¬ 
ments  in  the  form 


J  S«T  (X,t)  <r(X,t) 

Vo  I 


dV-  J  z  UT  (  X,»  )  F  (X,t)dV 
Vo  I 


-  J  8uT  (S,  t)  Fs(S,t)  dS  J  dt=  0 
Sr  ' 


(16) 
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where 


F  v(  X  ,t)  =  prescribed  body  force  vector 

F  s( $  *t)  =  prescribed  surface  tractions  on  portion  of  boundary,  ST. 

Su  (X  ,t)  =  [  S  v(x,y,t)  ]  =  x  and  y  directed  virtual  displacements 

S  denotes  the  boundary  position  vector. 


Linear  displacements  are  assumed  in  each  element.  This  leads  to  constant  strains  in  each 
element  giving  the  familiar  expression 
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cq  co 

xba 

1 

*< 
cr 
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«  (t  )  =  <j>  q  (t  ) 


ua(t) 

*o<t) 
Ub(t) 
vb(t) 
Uc(  t  ) 
Vc(t  ) 


(17) 


where  referring  to  Figure  3 


ybc  '  yb  -  V  *bc  ‘  *b  -  V  et0- 

A  =  element  area 


Linear  distribution  of  temperature  in  each  element  is  assumed,  giving 


®  (X,t  ) 


©o  (  t  ) 

©b(  t  ) 

©<:<♦> 


(18) 
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where  again  referring  to  Figure  3 


N 


2A 


xb*c  ~xc  *b 


-x 


be 

be 


c  jq  a  'c 


co 


xa  *b  "xb  ya 


-y 


ca 


ba 

;ba 


(  19) 


©a(t),  ©^(t),  ©c(t)  =  nodal  temperature  changes. 

The  shifted  time  is  calculated  using  an  average  nodal  temperature  in  each  element 

.  ^  4  Cx  h[Tm(t')ldt' 

i  --  t  (  t  )  =  J  10  L  m  J 


(20) 


where 


m 


(t)=(- 


T(X„  »  )+T  (Xh  t  )  +  T(X~  t  ) 


j1  average  nodal  temperature  for 
an  element. 


Insertion  of  Equations  14,  17,  and  18  into  Equation  16  results  in 
rz 

J  (§qT(U  (  /  L  c(X,t  )  +  C  q(t ) 

t,  Vo  I 

-3Kal  [!  xy]  N  0  (t  )  dV  -)  SqT(t  )  F(t  ))  dt  =  0 

where  the  virtual  work  of  body  and  external  loads  has  been  replaced  by  the  virtual  work  of  a 
statically  equivalent  set  of  nodal  forces.  Upon  integration  throughout  the  element  volume,  and 
assuming  the  Sq.(t)  are  independent, 


(21) 


where 


L  c  ( X  ,  t' )  +  <f?  C  <f>  q  ( t  ) 

-3Ka$T  X  [  I  xm  ym]  N  0  (t)  =  F(t  ) 

xa  +  xb  +xc  *o  +  *b  +  ^ 


(22) 


m 


ym 


3  3 

The  force  vector,  F  (t),  is  to  be  interpreted  as  the  forces  per  unit  element  thickness. 

The  product  3K  T I  jl  xm  ym  j  N  is  given  the  designation  D  .  Then 


3Ka 
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~xbc 
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-  *ca 
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(23) 
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where 

p  ixbyc  -xcyb  +xm  ybc  -  *m  xbc 

®  x  c  ya  -  xa  y c  xm  ^ca  -  x  co 

R  =  x„  yK  -  x  K  y„  —  x  m  y  K  „  F  y„,  x .  „ 

a  7  b  d  'a  m  ’  ba  7  m  Da 

T 

Letting  =  <j>  C  <f>  ,  the  resulting  expression 

<£T  L  t  (X,t)  F  K,q(t  )  -  D  01  t  )  =  Fit)  (24) 

looks  familiar  except  for  the  first  term,  which  requires  time  integration  from  t1  =  0  to  t'  ~  t. 
This  term  (see  the  one  dimensional  example  in  the  previous  section)  is  approximated 
numerically  as  follows  for  the  kth  time  point. 


<j>J  L  c  ( X,  t  )=  <f> 


T  I 


r-4  2  0  ] 

2-4  0 
0  0-3 


i=K-2 


(I  [3<«k-t|tl)-3<fk-0] 


+  ^-[g(o)-G 


>x(tk-D 
ey(tk  -I) 

Xy  {tk  -») 

€x  Hk  ) 

€y  (tk  ) 

l/./V 


rX(.i  )' 

e*  ( t  j  ) 

y*(tj  ) 
'xy  1 


(25) 


where  the  average  strain 


v  .  ~x 
€x  (t  i  ]  -  - 


It  j  +  j  )  +  ex  (tj) 


,  etc . 


All  of  the  strains  in  Equation  25  are  known  except  |^€x(tk)  yxy^)  j.  Using  Equation 

17,  we  write 


L  «  (X  ,t'  )  =  -  V  (  tk  )  F  [g(o  )  -  G  (£k  -Ck.  ,  )]  Kz  q  (tk) 


(26) 


where 


k  -2 


V(.k)  =  - m  (  X  [Gi€t-€i+,  >-G«fk  -e, )]«(»,  ) 

i  =  t 

+  <fk-fk-,']4 (,k-i  >) 


(27) 
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*2<»k 


and 


0 

0 

'3 


The  governing  incremental  stiffness  equations  for  an  element  are  now 


(k,  +[G(ol-G(£k-fk.,)]  K2),(tk) 
(Constant)  (Constant) 


=  F  (tk  )  -t-  H(tk)  +  V(t  k  ) 

(Meehan-  (Thermal  (Memory 
ical  Load  )  Load) 
Load  I 


(2  8) 


The  element  stresses  are  written  for  the  element  centroid  since  temperature  is  permitted  to 
vary  linearly  in  each  element. 


O’  x  (  t^  ) 

« x(tk  ] 

"®Q(tk  ) 

°*y  < f  k ) 

a  8,  +  S2 

€y  (tk  ) 

+  S3 

®b(tk> 

_rxy  (tk)  _ 

Yx  y(,k> 

W 

where 


s.  =  ? 


r-A  2  0 

2-4  0 

0  0-3 


K-2 


Vi 


[si«. 

I  =  I 

+  T  [G,°>  “G  <£k  >]  «<•»-, MVV 


-4  2  0 

2-4  0 

.  0  0-3 


3  K  q 
2  A 


P  Q 
P  0 
0  0 


R 

R 
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(29) 
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Using  the  development  in  the  preceding  section,  the  specialization  of  Equations  28  and  29  to 
Zak’s  method  is  accomplished  by 


k  ^k-i^ 


q 

1 

r  3-  r.  i 

V<V  ■  "  x  Ai  “]. 

j1' 

,  k  +  ~Z 

[g(o  )-G(cO)-£  a.  e  J  J 

jr' 

-(  £  k"£k 

-(  £k-i“£k-2^ 

°i.«  =  •  r| 

..  =  a .  - [ o  ] 

J  ,  '  J .  2  L  0  J 

[(' 

f. 

ri  )«*»<,-.  1  +  “i.» 

«i,3  -  ^  ' 

(•  r’ 

-  1 j  «  10) 

One  potential  advantage  of  Equation  28  should  be  noted.  Since  K  ^  and  K  ^  are  constant 
elemental  matrices,  for  isothermal  conditions  and  equal  time  increments,  the  merged  stiffness 
matrix  will  remain  constant  from  step  to  step.  And  for  isothermal  conditions  and  unequal  time 
steps,  a  merged  K  ^  and  merged  K  ^  can  be  stored  for  all  time,  thus  eliminating  regeneration 
and  merging  of  these  elemental  matrices  at  each  time  point.  Similar  advantages  may  be 
possible  for  transient,  homogeneous  temperature  distribution. 
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SECTION  IV 

THERMOVISCOELASTIC  ANALYSIS  PROGRAM  (TAP)  COMPUTER  CODE 


The  TAP  code  was  written  as  a  pilot  program  to  study  the  feasibility  of  our  approach. 
Originally,  out-of-core  storage  was  provided  for  so  that  an  unlimited  number  of  solution  time 
steps  could  be  used.  The  need  for  frequent  program  changes  coupled  with  excessive  turn¬ 
around  time  made  it  advisable  to  write  an  in-core  version  which  is  presently  the  operable 
program.  Coding  is  in  FORTRAN  IV  for  the  UNIVAC  1108  computer.  Zak’s  method  is  included 
as  an  option,  although  none  of  the  results  presented  in  the  figures  were  obtained  by  this  method. 

As  an  indication  of  the  program  capacity  and  efficiency,  the  problem  shown  in  Figure  10  for 
91  time  points  took  9  minutes  of  computer  time.  Approximately  5  minutes  of  this  time  was  spent 
in  generating  memory  loads  whereas  using  a  10-coefficient,  Prony  series  with  Zak’s  method, 
approximately  one  minute  is  necessary  for  generating  memory  loads.  For  this  problem, 
approximately  100  time  points  could  be  normally  used  before  exceeding  core  limitations  where¬ 
as  Zak’s  method,  the  number  of  time  points  is  essentially  unrestricted. 
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SECTION  V 

SOLUTION  OF  PROBLEMS 

ELEMENTARY  PROBLEMS 

In  order  to  evaluate  the  basic  operation  and  accuracy  of  the  TAP  program,  a  problem  with 
simple  geometry  was  chosen  for  which  exact  analytical  solutions  could  be  obtained.  The 
problem  geometry  is  illustrated  in  Figure  4.  The  material  properties  chosen  do  not  correspond 
to  any  particular  material. 


The  first  loading  (F-l)  considered  was  a  unit  step  load  of  10  on  the  right  hand  face.  The 
exact  solution  is  easily  found  to  be 


77  •  l0'S[ 


70.  3494  H  (  t  )  -  .0  4942  e 


-  2.0  588 1 


-  69.  3000  e 


-  .02398  t 


The  TAP  results  using  four  elements  are  shown  in  Figure  5  along  with  the  corresponding 
finite  element  solutions  by  Schapery’ s  direct  method  and  the  quasi- elastic  method.  The 
plotted  results  are  a  valid  comparison  between  the  three  methods  only  to  the  degree  that  the 
material  properties  represent  a  realistic  viscoelastic  material.  A  transform  parameter¬ 
time  correspondence  factor  of  p  =  -^-was  used  in  the  Schapery  method.  In  view  of  the  form  of 
the  relaxation  modulus  curve  in  Figure  1  the  quasi-elastic  result  was  not  unexpected. 


Next  considered  was  a  temperature  loading  (T-l), 

T  {x,t)  =  H(t)(l8  +  6x  )(  I  — .5  e~2t  }  t 

with  the  structural  and  material  reference  temperatures  taken  as  zero.  Following  essentially 
the  procedure  in  Reference  10  the  exact  solution  is: 


u  {  X  ,  t  } 


=  1.2857 


/"[Tlx'.l  )+  .  339le-|2-0588/,|0h[T1’'''  ’’  ’] 


f*  (2.0588  f  IOh[T(x'-t/')]dt")  hrT,  1 

Je  Jo  T<x',t')-IOhr(x*t  >Jdf' 

o 
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Figure  5.  Elementary  Problem,  F-l 
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Figure  6,  Elementary  Problem,  T-l 
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The  four  element  TAP  solution  (Figure  6),  using 0.1  time  intervals,  essentially  coincided  with 
the  exact  solution.  The  quasi-elastic,  finite  element  solution  was  nearly  exact  also  which  is 
not  surprising  in  view  of  the  different  type  of  loading  and  loading  rate  in  problem  T-2. 

COMPLEX  PROBLEMS 

Two  problems  were  chosen  representing  practical,  solid  rocket  grain  configurations, 
Figures  7  and  8.  The  elemental  breakdowns  are  obviously  not  fine  enough  for  a  realistic 
stress  analysis,  but  are  sufficient  for  evaluating  our  methods. 

The  first  problem  (F-2)  is  an  end  burning  grain,  partially  bonded  circumferentially. 
Figure  7  shows  a  one- sixth  segment  of  the  cross  section,  perpendicular  to  the  motor  axis. 
The  loading  represents  a  70°F,  firing  condition  for  which  a  step  pressure  and  step  radial 
displacement  is  applied  over  the  unbonded  and  bonded  portions  respectively.  Using  experimental 
data  in  Figure  8,  a  shear  modulus,  G(t)  =  1/3  E(t),  was  assumed.  Equal  increments  of  log 
time,  A  log  1Qt  =  1,  were  used.  Computer  running  time  was  approximately  2.5  minutes. 
Each  of  the  quasi-elastic  solutions  required  20  seconds.  The  TAP  results,  compared  with 
quasi-elastic  solutions  in  figure  7  demonstrate  the  inadequacy  of  the  quasi-elastic  method  for 
this  particular  problem. 

The  second  problem  was  a  slotted  grain  configuration.  Figure  10  shows  a  one-quarter 
segment  of  the  cross  section  perpendicular  to  the  motor  axis.  Assumptions  include  a  rigid 
outer  boundary,  a  uniform  structural  reference  temperature  T^  =  70°F,  and  a  shear  relaxation 
modulus  shown  in  Figure  2.  Tims-temperature  shift  data  from  Figure  9  determined  C  ±  and 
C  2-  A  temperature  (T-2),  representing  cooling  of  the  outer  surface  from  70°F  to  -70°F  was 
first  input  as  follows 

T(r,t)  =  70  -  70  (1  -  Cos  rr  /5t)(r/8)3  r  =  radial  location 

the  results  for  which  are  shown  in  Figure  10.  The  quasi-elastic  solution,  a  commonly  used 
method,  agrees  reasonably  well  for  displacements  but  deviates  greatly  from  the  TAP  solution 
for  stress.  The  TAP  results  appear  to  converge  with  progressively  finer  time  increments  as 
indicated  by  the  tabulated  solutions  in  Figure  10, 

A  slower  rate  of  cooling  was  considered  finally  as  follows: 

T(r,t)  =  70  -  70  (1  -  Cos  tt  /60t)(r/8)3 

It  was  anticipated  that  the  quasi-elastic  stress  solution  would  be  closer  to  the  TAP  results. 
Such  was  not  the  case  as  shown  in  Figure  11. 
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Figure  7.  Problem  F-2 
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SHIFT  FACTOR 


Figure  9.  Time-Temperature  Shift  Data  for  Problem  T-2, 
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Temperature  input:  T(r,t)  =  70  *  70(1-  Cos  t)(*g) 
Struct.  Ref.  Temp,  ;  T,  =  70* F 


TIME  ,  minutes 


Figure  11.  Slotted  Grain  Thermal  Problem  T-2A 
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SECTION  VI 
CONCLUSIONS 

The  two  primary  questions  under  study  are:  Within  the  assumptions  of  linear,  quasi¬ 
static  viscoelasticity,  what  is  the  adequacy  of  the  presently  used,  approximate  viscoelastic 
analysis  methods?  Is  the  direct  formulation  here  using  a  finite  element  procedure  a  practical 
approach  to  viscoelastic  stress  analysis? 

To  definitively  answer  the  first  question,  considerably  more  study  is  required.  In 
particular,  thermal  cyclingproblems  and  mechanical  loadingunder  thermal  environment  should 
be  studied  for  carefully  chosen  geometric  configurations.  Obviously,  the  question  must  be 
rephrased  if  real  materials  exhibit  significant  non-linearities,  as  appears  to  be  the  case  for 
composite  solid  propellants  for  example. 

The  answer  to  the  second  question  appears  to  be  in  the  affirmative.  In  regard  to  computer 
time  requirements,  problem  T-2  for  54  time  points  was  extrapolated  to  the  case  of  500  un¬ 
constrained  freedoms  (versus  60  in  T-2).  Based  on  the  TAP  code,  a  rough  estimate  of  2  hours 
computer  time  would  be  required — using  an  out-of-core  program  of  course.  Using  Zak’s 
method  in  an  in-core  program,  with  10  Prony  coefficients,  approximately  45  minutes  would 
be  required.  For  91  time  points  the  corresponding  figures  are  3  1/2  hours  and  with  Zak’s 
method,  80  minutes.  Significant  improvement  could  be  expected  with  a  more  efficient  code. 


For  short  solution  time  realms,  as  in  problem  F-2,  computer  time  would  not  generally 
be  prohibitive  for  realistic,  elemental  breakdowns.  Ofsignificanceinthiscaseis  the  practicality 
of  using  directly,  measured  relaxation  data  rather  than  resorting  to  additional  approximations. 
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